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Recent years have witnessed a rapid advancement of high-power short-pulse laser 
techniques which make it possible to observe many striking phenomena such as 
multi-photon ionization (MPI) and above-threshold ionization (ATI). This progress 
in experimental techniques has been accompanied by equally rapid development of 
theoretical methods needed to describe adequately phenomena occurring in strong laser 
fields. Representative reviews giving detailed account of these theoretical advances can 
be found in Protopapas et al (1997), Lambropoulos et al (1998), Chu & Telnov (2004) 
and Posthumus (2004). 

It is well-known that interaction of an atom and the electromagnetic (EM) field 
can be described in various ways corresponding to different choices of the gauge. In 
single-photon ionization calculations, it is the length and velocity gauges that are used 
most commonly. This choice is quite natural since in the length and velocity gauges the 
dipole matrix elements assume a very simple form. If, however, one wishes to compute 
probabilities of various multi-photon processes, a difficulty immediately arises. The 
dipole matrix elements corresponding to the free-free electron transitions are divergent. 

There are various ways to circumvent this difficulty For one-electron systems 
this problem can be easily avoided since analytical expressions for the Coulomb Green 
function are known. This fact has been exploited in many papers (Klarsfeld 1969, 
Rapoport et al 1969, Arnous et al 1973, Teague & Lambropoulos 1976, Klarsfeld & 
Maquet 1979, Karule 1971, Karule 1985, Karule & Moine 2003). Alternatively, one 
may avoid the divergency problems by reducing summations and integrations over 
intermediate electron states to solution of nonhomogeneous differential equations (Zernik 
& Klopfenstein 1965, Chan & Tang 1969, Jayadevan & Thayyullathil 2001). 

For systems with more than one electron, where these techniques cannot be 
implemented, other methods have been developed. The states belonging to the 
continuous spectrum of the system can be represented by means of a suitable set 
of the square integrable (L 2 ) functions (Tang & Bachau 1993, Venuti et al 1996). 
The system is quantized in a box of sufficiently large dimensions which gives a 
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discretized representation of the continuous spectrum. This approach, supplemented 
by the B-splines technique, allowed to obtain a set of accurate MPI and ATI 
cross-sections for many-electron systems such as He (Venuti et al 1996, Saenz & 
Lambropoulos 1999, Nikolopoulos & Lambropoulos 2001) and Be (Tang & Bachau 1993). 
Another technique which was applied successfully to study MPI in many-electron 
systems uses a regularization procedure for the free-free matrix-elements (Mercouris 
& Nicolaides 1989, Nicolaides & Mercouris 1989, Korol 1997). 

In the present paper we describe an alternative method allowing to compute 
amplitudes of the MPI and ATI processes in many-electron systems in a more direct way. 
The proposed method is based on the so-called Kramers-Henneberger (KH) form of the 
interaction Hamiltonian (Kramers 1956, Henneberger 1968). The KH representation 
of the system "atom plus electromagnetic field" is often used when one studies a 
temporal evolution of atomic systems subjected to a pulse of EM radiation (Reed & 
Burnett 1991, Vivirito & Knight 1995). The KH representation also enables simple 
asymptotic boundary conditions used in the external region of the i?-matrix Floquet 
theory (Burke et al 1991). 

In the present paper we shall be interested in another aspect of the KH 
representation, namely the advantages its use may offer in the perturbative 
computations of MPI rates. To our knowledge, the KH description of the EM radiation 
interaction with atoms has not been used in this context. As we shall see, in the 
perturbative calculations the KH representation offers one important advantage. In 
contrast to the length and velocity gauges, in the KH formulation all the dipole 
matrix elements are finite and well-defined. For the laser fields of not very large 
intensities (below 10 13 Wcm" 2 ) the perturbation theory (PT) provides quite an adequate 
description of the MPI process and allows to achieve numerically accurate results with 
much less computational labor. We shall consider below two-photon ionization processes 
in the hydrogen and helium atoms. The highly accurate perturbative results available in 
the literature allow us to evaluate directly the accuracy of the method. We also discuss 
some subtle numerical aspects of application of the KH Hamiltonian in perturbative 
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calculations. 

The use of the KH Hamiltonian in time-independent MPI calcualtions is not 
restricted to preturbative regime. In our recent paper (Ivanov & Kheifets 2005) we 
demonstrated utility of this method in a non-perturbative calculation of the MPI rates 
in atomic hydrogen. 

In the following sections, we briefly recall main theoretical aspects of the KH 
representation, describe the computation of the dipole matrix elements and apply the 
technique to two-photon ionization of atomic hydrogen and helium. 

2. Theory 

A starting point for the description of the interaction of an atom and a monochromatic 
EM field is the minimal coupling Hamiltonian: 

H = i?atom + H int . (1) 

Here H a t m has the usual meaning of the atomic Hamiltonian: 

N p 2 N Z N 1 

#atom = ^lf-^-+ ]T — (2) 

i=\ * i=l ' * i,j=l,i>j 1 V 

The atomic Hamiltonian is taken in a non-relativistic form with Z being the nucleus 
charge. The part of the Hamiltonian H int which describes the interaction of the atom 
and the EM field can be written as (see e.g. Sobelman (1972)) 

i N ( ~ aA 

ffint = __g^.p,__j, (3) 

where p is the momentum operator, A is a vector potential, summation runs over all 

atomic electrons. In the following we shall assume that the dipole approximation is 

valid so that the vector potential does not depend on atomic coordinates. Performing a 

suitable canonical transformation of Equation (3), one can obtain various forms of the 

interaction Hamiltonian. The KH Hamilton is obtained by the canonical transformation 

-^kh = e 1 ^ H min e~ 1 ^ — generated by the operator: 

at 

-i * N 1 * 

f = — J £ A{r) Pl dr+—J A\r) dr , (4) 

n 1=1 n 
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Expressed in quantum-mechanical terms (as far as description of the EM field 
is concerned), this transformation is also known as the Pauli-Fierz canonical 
transformation (Pauli & Fierz 1938). We shall not distinguish between these two versions 
of the transformation as the final results are identical. 

Under the transformation (4) the minimal-coupling Hamiltonian (2) becomes: 

HkH = -f^atom + #™ ? (5) 

where H&tom retains the same form as the Hamiltonian (2) while the interaction 
Hamiltonian becomes: 

AS" = E (t - tA^\ • (6) 



i=l 



Ti Vi + OL 



We shall be interested in the case of a linearly polarized monochromatic EM field. In 
this case a = F/lu 2 where F is the operator of the electric field intensity, u is the photon 
energy. §. If a can be considered as a small quantity, the leading term of expansion 
of (6) reproduces the well-known Zr/r 3 form of the interaction Hamiltonian in the 
acceleration gauge. This form is often used in the first-order perturbation calculations. 
We, however, are interested in higher order effects and must generally retain complete 
form of the Hamiltonian (6). 

To build the perturbation theory expansion, treating operator (6) as a perturbation, 
we need a formula for the matrix elements of this operator sandwiched between the 
states describing the noninteracting atom and the EM field. It is convenient for our 
purposes to use the notation \a, m) for these states where a stands for a set of quantum 
numbers describing the atom and m denotes a number of laser photons in a given mode. 
Such a formula can be obtained from the known matrix elements of the quantized 
vector potential operator (Ivanov & Kheifets 2005). A simpler derivation, relying on 
the correspondence between the quantum and classical description of the EM field, is 
given in the Appendix I. 

Obtained in either way, the final formula reads: 

' Z Z 



rVKH 



^a,n + p 

§ We use the atomic system of units in which h = e = m 



1 N J I 
b, nj = — ^2 / cos p9 / a 

71 i=l o \ 



\ri + Fcos#/u; 2 | 



b ) d0(7) 
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Here F is already a classical vector and not an operator. Its magnitude is related to 
the number of photons via the flux conservation relation F 2 /8tt = ncu, and it is directed 
along the polarization vector of the incoming photons. In the following, we shall take 
this direction as the z-axis. 

To amend Equation (7) to a form suitable for practical computations, we use the 
well-known expansion: 



where r< (r>) is the smaller (greater) of r and F cos 9/uj 2 . Equation (8) allows separation 
of the radial and angular variables. Angular parts are evaluated analytically using 
integrals of products of several spherical functions (Varshalovich et al 1988): 



Equation (9) is written for the case when the atomic subsystem contains one electron 
(hydrogen). Generalization for the case of many-electron systems is a simple exercise in 
angular momentum algebra. 

When performing perturbation calculations, we are usually interested in keeping 
track of field dependencies of the matrix elements. Suppose we study a process for 
which, in the leading order, the amplitude is proportional to k-th power of electric 
field strength M oc F k , with some integer k. Then we would like to retain in the 
perturbation theory expressions only the terms which give rise to such dependence in 
the limit of small F . If the length or velocity forms for the atom-EM field interaction 
is used, such a count of powers of F is trivial, following from the well-known selection 
rules for the matrix elements. In the KH representation, there are no exact selection 
rules. Nevertheless, the count of powers of the electric field strength is still possible. 
Consider, for example, the case of a two-photon ionization of a one-electron atom, which 
is the second-order process. Suppose, we are interested in the ionization from the state 
a with a given orbital momentum /. Then it is easy to see from Equation (8) that the 
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following asymptotics holds for F — > 0: 

(al m\H^\bl m ± 2) oc F 2 , (al m\Hgj*\bl ± 1 m ± 1) oc F . (10) 

Here a and b stand for the set of all atomic quantum numbers except the angular 
momentum, the integer m refers to a total number of laser photons. It is easy to see 
from Equation (8) that the coefficient of proportionality in the second matrix element 
of (10) is just a matrix element of the operator Zr/r 3 which is commonly used in the 
first order calculations in the acceleration gauge. 

Thus, the leading terms of the amplitude of the two-photon ionization can be 
written as: 

M(al m^bl m-2) = (al m\H^\bl m - 2) 

i — l)(cl' m — 1 
E a i — E d > + ie 



+ v£ (oj m \ H ™\cl' m ~ l)(c*' m - l|#jnt H |^ m-2) ^ 



cV l'=l±l 

M(al m -> bl + 2 m - 2) = (al m\H™\bl + 2 m-2) 

^ (al m\H™\cl' m-l)(cl' m - l\H™\bl + 2 m - 2) 

+ „ h , ^ - ^ + * e 

d' i'=Z+l 

M(aZ m ^ bl -2 m-2) = (al m\H™\bl -2 m-2) 

y-jk (al m\H^\cl' m- l)(cl' m- -2 m-2) 

Here the symbol ^ indicates the sum over the discrete spectrum and integration over 
continuum of the intermediate states. Without sacrifice of accuracy, in these sums we 
may use the operator Zr/r 3 instead of the complete form of the operator H™. 

Once the amplitudes (11) - (13) are computed, the generalized partial cross-section 
of the two-photon ionization from the initial state a, I to a final channel b, I' is given by 
(c.f. Tang & Bachau (1993)) 

a(al - bl') = 2 7 n 3 a 2 a 4 r u 2 lim \M^Lz^DL [cm 4 s -i] (14) 

f^o F 4 k 

Here a, do, To are the fine structure constant, the atomic unit length in cm and the 
atomic unit time in seconds. F is the EM field strength and oo is the photon energy, 
both expressed in the atomic units. The one-electron continuum wave functions used 
to calculate the ionization amplitudes are normalized on the momentum scale. The 
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generalized cross-section (14) is related to the ionization rate 

r(ai - bi>) = a{al hV) x 1Ql2A Ryd [W" W] , (15) 

V ; 13.605 x 1.60219 L J ' y ! 

where A is wavelength (in nm) and Ryd = 109677 cm -1 is the Rydberg constant (Karule 
& Moine 2003). 

It has to be noted that the KH transformation modifies, in general, the atomic states 
after the field is switched off at a sufficiently large time (Vivirito & Knight 1995, Reed & 
Burnett 1991). However, it is shown in Appendix II that this does not affect the cross- 
sections or ionization rates and can only be noticed in the fine details of the temporal 
evolution of the system 

3. Numerical computations 

Certain amount of care has to be exercised when the amplitudes (11) - (13) are computed 
numerically. Consider matrix elements in the sums over intermediate states in these 
expressions. Separating radial and angular variables with the help of equations (7), 
(8), and taking the form Zr/r 3 for the interaction Hamiltonian we obtain integrals of 
the sort I = J REi(r)Rki'(r) dr. Here -R's are radial electron wave functions behaving 
as r l near the origin r = 0. Function R Ei describes either an initial or final state of 
the process, Ruv is the radial wave function of the intermediate state belonging to the 
continuous spectrum which, for large momenta, behaves as Rkv(j : ) oc sin(A;r + S)/r. 
The /c-dependence of the integral is crucially dependent on the orbital momentum I. 
If / > the integral I can be approximated for large k as / Rei{t) sm.{kr + S)/r (we 
omit unimportant normalization factors). Consider the function P(r) = Rei/t. If / > 0, 

oo 

/ \P(r)\dr is finite, hence the Riemann-Lebesgues lemma is applicable and the integral / 

o 

will decay for k — > oo. If I = 0, the replacement of the Ruv by its asymptotic expression 
is not legitimate (we would obtain a divergent integral). More careful study shows that if 
/ = than / — > const when k — > oo. The integrals over momenta of intermediate states 
in the formulas (11) - (13) still converge in this case due to the energy denominator but 
much more slowly than in the case of / > 0. This means that to achieve a good numerical 
accuracy for the amplitudes with 1 = 0, one must take into account an asymptotic tail 
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of the integrand for k — > oo. 

In practice, this does not pose serious difficulties. One has only to determine (either 
analytically or numerically) the constant in the expression / — > const for k large enough 
so that this asymptotic law holds well, and then add the corresponding contribution to 
the integrals in formulas (11) - (13). For / > one need not worry about the asymptotic 
tails of the integrand as integrals / decay quite fast with k. 

For the presently considered targets (H and He), the initial state has an S orbital 
character leading to the two final channels S and D which correspond to amplitudes 

(11) and (12). As explained above, we may retain in these amplitudes only the P 
intermediate states. Both the discrete and continuous intermediate states were taken 
into account. The continuous spectrum integration follows closely prescriptions given in 
Bray (1994) and Ivanov & Kheifets (2005). The interval of momenta (0, g max ) is divided 
into several subintervals. For the photon energies above the ionization threshold, a 
pole is present in the momentum integral. To carry out the integration around the 
pole accurately, the first two subintervals are chosen to be (0, fc po ie) and (k po \ e , 2k po \ e ) 
with a typical number of 20 momentum points in each subinterval. Then the delta- 
function singularity is isolated and the remaining principle value integral is evaluated 
by a modified Gaussian rule (Bray 1994). The remaining part of the momentum integral 
is divided as follows: (2k po i e , 4) (20 integration points), (4, 10) (20 points) and (10, g max ) 
(20 points). These intervals are pole-free and the integration is performed by using a 
Gauss quadrature rule. The fairly large value of g max is chosen to take care of a slow 
decay of the integrand in the S'-channel. The asymptotic tail (g max , oo) is calculated 
analytically. 

It is worth to be noted that the first order matrix elements in the amplitudes (11), 

(12) are roughly of the same magnitude as the second order terms and, sometime, of 
the opposite sign. So their inclusion is essential. 

Once the amplitudes are computed the partial cross-sections and rates can be 
determined via Eqs. (14), (15). Total cross-sections and rates corresponding to linear 
and circular polarization of the EM field can then be determined. If inital state of atomic 
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system is an 5-state, than: T l (ns) = T(ns — > ks) + T(ns — > kd) (linear polarization) 
3 

and T c (ns) = -T(ns — > fed) (circular polarization) (Jayadevan & Thayyullathil 2001). 

All calculations reported below were performed for the EM field strength of 
F = 0.03 a.u. This field strength is small enough, so that determination of the cross 
sections via Eq.(14) for small but finite F will be, as we shall see, quite accurate. 

4. Results 

4-1. Hydrogen 

In the case of hydrogen, all the one-electron wave functions are know analytically. These 
allowed us to include 15 discrete intermediate P states and carry out integration up to 
Qmax = 70 a.u. in the second-order amplitudes. 

Table 1. Ionization rates (in units of W _1 cm 4 ) for the two-photon ionization of 
atomic hydrogen in the ground state by linearly T l and circularly T c polarized light. 
Numbers in parenthesis indicate powers of 10. The literature values are denoted as 
JT01 (Jayadevan & Thayyullathil 2001) and KM03 (Karule & Moine 2003) 



A 










r c 




nm 


JT01 


KM03 


Present 


JT01 


KM03 


Present 


20 


3.02(-38) 


3.01(-38) 


3.02(-38) 


2.44(-38) 


2.43(-38) 


2.44(-38) 


40 


2.15(-36) 


2.14(-36) 


2.15(-36) 


2.03(-36) 


2.02(-36) 


2.03(-36) 


60 


2.62(-35) 


2.61(-35) 


2.62(-35) 


2.79(-35) 


2.78(-35) 


2.79(-35) 


80 


1.58(-34) 


1.57(-34) 


1.59(-34) 


1.85(-34) 


1.84(-34) 


1.86(-34) 



In Table 1, we present the rates of the two-photon ionization of the ground state 
hydrogen atom with the linear and circular polarized light and compare them with the 
latest literature values (Jayadevan & Thayyullathil 2001, Karule & Moine 2003). The 
three sets of data are virtually identical for all the photon energies considered. The 
cited literature values have been obtained analytically which is possible for hydrogen. 
They can be, therefore, considered as "exact". The comparison in the Table indicates 
that the use of the KH Hamiltonian allowed us to achieve comparable level of accuracy 
with little effort. 
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4-2. Helium 

We consider now two-photon ionization from the ground state of helium. Due to possible 
two-electron excitations, complete treatment of this problem is much more complicated 
than for hydrogen. However, for the photon energies below the threshold of the N = 2 
excitations, we can use a frozen core approximation and restrict ourselves with only one 
active electron. In this approximation, the problem is effectively reduced to a hydrogen- 
like calculation with numerical Hartree-Fock wave functions. 

We describe the helium atom as follows. For the ground state we use the self- 
consistent Hartree-Fock approximation (Chernysheva et al 1976). The ground state is 
thus represented as a product of the two Is orbitals. For the excited states, both discrete 
and continuous, we use the frozen core Hartree-Fock approximation (Chernysheva 
et al 1979). These states are thus represented as properly symmetrized products of 
the Is core orbital and an orbital describing excited electron either in the bound state 
or the continuum. 

Because of this frozen-core approximation, the calculation for helium proceeds in 
almost exactly the same way as for hydrogen. The only difference is that due to the 
equivalence of the Is 2 electrons an additional factor of 2 arises in the formula for the 
cross-section. As for hydrogen, we may retain only P states in the sums over the 
intermediate states. We have retained 7 P-states corresponding to the excitations lsnp 
with n = 2 — 8. Integration over the continuous spectrum was performed as for hydrogen 
except for the value of g max = 30 a.u. All we said above concerning importance of the 
correct account of the asymptotic tail in the momentum integral for the S'-wave, applies 
for the case of helium as well. 

Our results are presented in the Figure in comparison with other calculations (Saenz 
& Lambropoulos 1999, Nikolopoulos & Lambropoulos 2001) which use a considerably 
more accurate representation of the helium atom. Despite a rather crude character 
of the frozen-core approximation, in the region of the photon energies considered we 
achieve quite a satisfactory agreement with the literature values. We could not, of 
course, extend our calculation into the region of larger photon energies since there the 
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D channel S channel 




Photon energy CO, eV Photon energy CO, eV 



Figure 1. Cross section of the two-photon ionization from the ground state 
of helium Comparison is made with literature values marked as SL99 (Saenz & 
Lambropoulos 1999) and NL01 (Nikolopoulos & Lambropoulos 2001) 



5. Conclusion 

We have shown that the KH description of the atom-EM radiation interaction can be 
used efficiently in calculations of MPI and ATI processes in realistic atomic systems. 
The fact that the dipole matrix elements between continuous electron states are finite 
and well-defined quantities makes the calculation relatively simple, both numerically 
and conceptually. 

For hydrogen, the present results agree completely with those obtained in Jayadevan 
& Thayyullathil (2001) and Karule & Moine (2003). As the latter results are analytical 
and can therefore be considered as virtually exact, we may be confident that the present 
approach allows to achieve quite a high accuracy. Our method can also be applied to 
the systems with more that one electron, as the comparison with the data for the 
two-photon ionization of helium shows. Despite the fact that we used rather a crude 
description of the field-free helium atom (we omitted the core-excitation effects), we 
obtained good quantitative agreement with the results of other authors who employed 
a more elaborate representation for the helium atom. As was indicated above, the 
accuracy of our description of MPI of helium can be further improved without any 
problems of conceptual character. All we have to do is to "thaw" the core and to 
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allow the two-electron excitations. This can be done, for example, with the use of the 
convergent close-coupling (CCC) method (Bray 1994) which is known to provide good 
description of a complete set of two-electron states, both discrete and continuous. Such 
calculation will be reported elsewhere. 

6. Appendix I 

We give below a derivation of Equation (7) for the matrix elements of the KH interaction 
Hamiltonian based on the well-known correspondence between the classical Floquet and 
the quantum-mechanical descriptions of the atom-EM field interaction (Shirley 1965). 

In the classical picture, the KH Hamiltonian has the form (6) with a = F coscut/co 2 
where F is a classical amplitude of the EM field. With this expression being a periodic 
function of time, the Shcrodinger equation has a set of solutions (the Floquet anzats) 
which allows the following Floquet-Fourier expansion: 

+oo 

m(t) = e ~ lEt u ne muJ \ (16) 

where E is the quasi-energy. Expanding the time-periodic function as a Fourier 
series: 



H™(t) = Y,V n cosnut, (17) 

n=0 



with 



V n = — / H™(wt) cosnujt dt, (18) 

7T J 


and equating coefficients with e~' imuJt , one obtains a set of equations for the Fourier 
amplitudes u n : 

(E - nU) - ffatom) U n = 7«m + T^m (19) 

k,m,k>0 ^ k,m,k>0 ^ 

m—k—n m+k=n 

In the quantum-mechanical description, the coefficients with the amplitudes u m on the 
r.h.s of Equation (19) are nothing but the matrix elements (n\H^\m) taken between 
the states with n and m photons (Shirley 1965). This correspondence holds if we neglect 
all spontaneous processes and retain only laser photons. Since in Eq.(19) summation 
index k > 0, for given n,m, n ^ m, the two terms on the r.h.s of this equation can 



Kramers- Henneberger calculations of multi-photon ionization 



14 



be combined to give rise to 




This gives immediately the formula (7) for the 



2 



matrix elements of the operator ff™. In our earlier paper (Ivanov & Kheifets 2005), 
we obtained this formula directly using quantized form of the electric field operator. 

7. Appendix II 

A comment has to be made on the applicability in the present case of the Fermi golden 
rule which is used to derive Eq.(14). Below, all the discussion uses only classical terms 
for the description of the EM field which leads to the same results as the full quantum- 
mechanical treatment but is somewhat simpler. 

Physically, transformation to the KH frame is equivalent to transformation to the 
non-inertial frame oscillating with the electron. Let us suppose that an initially field- 
free atom is in some state \l/o- At the moment t = the interaction of the atom and 
EM field is switched on and it is switched off at the moment t — t 1 with some ti large 
enough so that all the transient processes are negligible. With the help of formulas (11) - 
(13) we can obtain (for sufficiently large ti) perturbative solution of the time-dependent 
Schrodinger equation in the KH- frame. As a result, we obtain at the moment t\ a vector 
in the KH frame: 



where ?7 KH (0,ti) is the evolution operator which we constructed in the KH frame to a 
given order of the perturbation theory. 

To determine rates of various processes, we must now find overlaps of the vector thus 
obtained with various states of the field-free atom. The latter, however, are generally 
different from the eigenstates of ifo, the original field-free atomic Hamiltonian. As 
formula (4) shows, they are connected to the eigenstates of if by means of a unitary 
transformation generated by the operator T in Eq. (4). The fact that final states are 
generally altered by this operator may play an important role in calculations of time- 
evolution of an atomic system under the action of an impulse of EM radiation (Vivirito 
& Knight 1995, Reed & Burnett 1991). 

The problem we are considering here differs from the problems of determining the 



*kh(*i) = U KH (0,i0*, 



0- 



(20) 
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time-evolution under the action of an impulse of EM radiation in one important aspect. 
We are at liberty to switch off interaction at any time (provided it is large enough). 
Operator T in Eq.(4) contains two terms. The one proportional to the square of the 
vector potential is unimportant (it is a pure phase factor in the dipole approximation 
which we use). It is the term linear in A that is responsible for the alteration of the 
final states. We may use, however, periodicity of the vector potential and choose the 
moment t\ so that this term disappears. The final states then remain unaltered and the 
validity of the Fermi golden rule is restored. 
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